{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Comparison Of Various Code Optimization Methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are a number of ways to optimize the performance of Python code.  Below is a short summary adapted from <a>http://people.duke.edu/~ccc14/sta-663/MakingCodeFast.html</a> regarding various optimization strategies.\n",
    "\n",
    "There is a traditional sequence for writing code, and it goes like this:\n",
    "\n",
    "1) Make it run<br />\n",
    "2) Make it right (testing)<br />\n",
    "3) Make it fast (optimization)<br />\n",
    "\n",
    "Making it fast is the last step, and you should only optimize when it is necessary.  Also, it is good to know when a program is “fast enough” for your needs.  Optimization has a price:\n",
    "\n",
    "1) Cost in programmer time<br />\n",
    "2) Optimized code is often more complex<br />\n",
    "3) Optimized code is often less generic<br />\n",
    "\n",
    "However, having fast code is often necessary for statistical computing, so we will spend some time learning how to make code run faster.  To do so, we need to understand why our code is slow.  Code can be slow because of differnet resource limitations:\n",
    "\n",
    "CPU-bound - CPU is working flat out<br />\n",
    "Memory-bound - Out of RAM - swapping to hard disk<br />\n",
    "IO-bound - Lots of data transfer to and from hard disk<br />\n",
    "Network-bound - CPU is waiting for data to come over network or from memory (“starvation”)<br />\n",
    "\n",
    "Different bottlenekcs may require different approaches.  However, there is a natural order to making code fast:\n",
    "\n",
    "1) Cheat\n",
    "* Use a better machine (e.g. if RAM is limiting, buy more RAM)\n",
    "* Solve a simpler problem (e.g. will a subsample of the data suffice?)\n",
    "* Solve a different problem\n",
    "\n",
    "2) Find out what is slowing down the code (profiling)\n",
    "* Using timeit\n",
    "* Using cProfile\n",
    "* Using memory_profiler\n",
    "\n",
    "3) Use better algorithms and data structures\n",
    "\n",
    "4) Off-load heavy computations to numpy/scipy\n",
    "\n",
    "5) Use compiled code written in another language\n",
    "* Calling code written in C (ctypes, cython)\n",
    "* Calling code written in Fotran (f2py)\n",
    "* Calling code written in Julia (pyjulia)\n",
    "\n",
    "6) Convert Python code to compiled code\n",
    "* Using numexpr\n",
    "* Using numba\n",
    "* Using cython\n",
    "\n",
    "7) Write parallel programs\n",
    "* Ahmdahl and Gustafsson’s laws\n",
    "* Embarassinlgy parallel problems\n",
    "* Problems requiring communication and syncrhonization\n",
    "\n",
    "8) Execute in parallel\n",
    "* On multi-core machines\n",
    "* On multiple machines\n",
    "* On GPUs\n",
    "\n",
    "This notebook will focus on 4 and 6.  We will use the example of calculating the pairwsise Euclidean distance between all points.  Examples are adapted from <a>http://people.duke.edu/~ccc14/sta-663/Optimization_Bakeoff.html</a>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "%precision 2\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import numexpr as ne\n",
    "from numba import jit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1000, 3)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xs = np.random.random((1000, 3))\n",
    "xs.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the pure python version of the algorithm.  Used as a baseline for comparison."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def pdist_python(xs):\n",
    "    n, p = xs.shape\n",
    "    D = np.empty((n, n), np.float)\n",
    "    for i in range(n):\n",
    "        for j in range(n):\n",
    "            s = 0.0\n",
    "            for k in range(p):\n",
    "                tmp = xs[i,k] - xs[j,k]\n",
    "                s += tmp * tmp\n",
    "            D[i, j] = s**0.5\n",
    "    return D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10 loops, best of 3: 3.17 s per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit -n 10 pdist_python(xs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numpy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NumPy is the fundamental package for scientific computing with Python. It contains among other things:\n",
    "\n",
    "* a powerful N-dimensional array object\n",
    "* sophisticated (broadcasting) functions\n",
    "* tools for integrating C/C++ and Fortran code\n",
    "* useful linear algebra, Fourier transform, and random number capabilities\n",
    "\n",
    "Besides its obvious scientific uses, NumPy can also be used as an efficient multi-dimensional container of generic data. Arbitrary data-types can be defined. This allows NumPy to seamlessly and speedily integrate with a wide variety of databases.\n",
    "\n",
    "Library documentation: <a>http://www.numpy.org/</a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def pdist_numpy(xs):\n",
    "    return np.sqrt(((xs[:,None,:] - xs)**2).sum(-1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100 loops, best of 3: 50.3 ms per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit -n 100 pdist_numpy(xs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numexpr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Numexpr is a fast numerical expression evaluator for NumPy. With it, expressions that operate on arrays (like \"3\\*a+4\\*b\") are accelerated and use less memory than doing the same calculation in Python.\n",
    "\n",
    "In addition, its multi-threaded capabilities can make use of all your cores, which may accelerate computations, most specially if they are not memory-bounded (e.g. those using transcendental functions).\n",
    "\n",
    "Library documentation: <a>https://github.com/pydata/numexpr</a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def pdist_numexpr(xs):\n",
    "    a = xs[:, np.newaxis, :]\n",
    "    return np.sqrt(ne.evaluate('sum((a-xs)**2, axis=2)'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100 loops, best of 3: 19.4 ms per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit -n 100 pdist_numexpr(xs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numba"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Numba gives you the power to speed up your applications with high performance functions written directly in Python. With a few annotations, array-oriented and math-heavy Python code can be just-in-time compiled to native machine instructions, similar in performance to C, C++ and Fortran, without having to switch languages or Python interpreters.\n",
    "\n",
    "Numba works by generating optimized machine code using the LLVM compiler infrastructure at import time, runtime, or statically (using the included pycc tool). Numba supports compilation of Python to run on either CPU or GPU hardware, and is designed to integrate with the Python scientific software stack.\n",
    "\n",
    "Library documentation: <a>http://numba.pydata.org/</a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "pdist_numba = jit(pdist_python)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100 loops, best of 3: 11.1 ms per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit -n 100 pdist_numba(xs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cython"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cython is an optimising static compiler for both the Python programming language and the extended Cython programming language. It makes writing C extensions for Python as easy as Python itself.\n",
    "Cython gives you the combined power of Python and C to let you:\n",
    "\n",
    "* write Python code that calls back and forth from and to C or C++ code natively at any point\n",
    "* easily tune readable Python code into plain C performance by adding static type declarations\n",
    "* use combined source code level debugging to find bugs in your Python, Cython and C code\n",
    "* interact efficiently with large data sets, e.g. using multi-dimensional NumPy arrays\n",
    "* quickly build your applications within the large, mature and widely used CPython ecosystem\n",
    "* integrate natively with existing code and data from legacy, low-level or high-performance libraries and applications\n",
    "\n",
    "Library details here: <a>http://cython.org/</a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%load_ext Cython"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%cython\n",
    "\n",
    "import numpy as np\n",
    "cimport cython\n",
    "from libc.math cimport sqrt\n",
    "\n",
    "@cython.boundscheck(False)\n",
    "@cython.wraparound(False)\n",
    "def pdist_cython(double[:, ::1] xs):\n",
    "    cdef int n = xs.shape[0]\n",
    "    cdef int p = xs.shape[1]\n",
    "    cdef double tmp, d\n",
    "    cdef double[:, ::1] D = np.empty((n, n), dtype=np.float)\n",
    "    for i in range(n):\n",
    "        for j in range(n):\n",
    "            d = 0.0\n",
    "            for k in range(p):\n",
    "                tmp = xs[i, k] - xs[j, k]\n",
    "                d += tmp * tmp\n",
    "            D[i, j] = sqrt(d)\n",
    "    return np.asarray(D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100 loops, best of 3: 11 ms per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit -n 100 pdist_cython(xs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Scipy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Scipy has an optimized version of this particular function already built in.  It exploits symmetry in the problem that we're not taking advantage of it in the \"naive\" implementations above.\n",
    "\n",
    "Library documentation: <a>http://www.scipy.org/</a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy.spatial.distance import pdist as pdist_scipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100 loops, best of 3: 5.79 ms per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit -n 100 pdist_scipy(xs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's all of them together."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Python\n",
      "10 loops, best of 3: 3.15 s per loop\n",
      "Numpy\n",
      "100 loops, best of 3: 49.5 ms per loop\n",
      "Numexpr\n",
      "100 loops, best of 3: 19.3 ms per loop\n",
      "Numba\n",
      "100 loops, best of 3: 11 ms per loop\n",
      "Cython\n",
      "100 loops, best of 3: 11 ms per loop\n",
      "Scipy\n",
      "100 loops, best of 3: 5.68 ms per loop\n"
     ]
    }
   ],
   "source": [
    "print('Python')\n",
    "%timeit -n 10 pdist_python(xs)\n",
    "print('Numpy')\n",
    "%timeit -n 100 pdist_numpy(xs)\n",
    "print('Numexpr')\n",
    "%timeit -n 100 pdist_numexpr(xs)\n",
    "print('Numba')\n",
    "%timeit -n 100 pdist_numba(xs)\n",
    "print('Cython')\n",
    "%timeit -n 100 pdist_cython(xs)\n",
    "print('Scipy')\n",
    "%timeit -n 100 pdist_scipy(xs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some observations:\n",
    "\n",
    "* Pure python is much, much slower than all of the other methods (close to 1000x difference!)\n",
    "* Simply using Numpy where possible results in a huge speed-up\n",
    "* Numba is surprisingly effective given how easy it is to utilize, on par with compiled C code using Cython\n",
    "* Algorithm optimizations (such as those employed in the Scipy implementation) can easily trump other methods"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
